## Plot daily average heat rate
## These results are reported in Figure 7
rm(list=ls())
options(scipen=100,digits=10)

library(dplyr)
library(ggplot2)
# Load data
load("data_final.RData")
Plant_List=unique(data3$Plant_Code)

cems1 <- read.csv("cems_112913_022714.csv")
cems2 <- read.csv("cems_112914_022715.csv")
cems3 <- read.csv("cems_112915_022716.csv")
cems4 <- read.csv("cems_112916_022717.csv")
cems5 <- read.csv("cems_112917_022718.csv")
cems6 <- read.csv("cems_112918_022719.csv")

cems7 <- read.csv("cems_112919_022720.csv")
cems8 <- read.csv("cems_112920_022721.csv")

cems <- rbind(cems1, cems2, cems3, cems4, cems5, cems6)

data <- cems %>%
  filter(Facility.ID %in% Plant_List) %>%
  group_by(Facility.ID) %>%
  mutate(heat.rate = Heat.Input..mmBtu. / Gross.Load..MWh.)

data <- data %>%
  group_by(Facility.ID) %>%
  mutate(units = length(unique(Unit.ID))) %>%
  filter(units > 1)

# create a dataset that includes Facility.ID and number of units
units <- data %>%
  group_by(Facility.ID) %>%
  distinct(units) 

data$numeric_date <- as.numeric(gsub("-", "", data$Date))

# calculate daily average heat rate
hr_avg <- data %>%
  filter(!is.infinite(heat.rate)) %>%
  group_by(numeric_date) %>%
  summarise(avg.hr = mean(heat.rate, na.rm = TRUE))


# convert numeric date to character date
hr_avg$date <- as.Date(as.character(hr_avg$numeric_date), format = "%Y%m%d")


# Set the starting day for the tick labels. Adjust dates to produce heat rate trend for different years.
start_day <- as.Date("2018-11-29")
# Generate a sequence of dates spaced by 7 days starting from start_day
tick_dates <- seq(from = start_day, length.out = 14, by = "7 days")
# plot daily average heat rate 
ggplot(hr_avg, aes(x = date, y = avg.hr)) +
  geom_line() +
  labs(x="Date", y = "Average Daily Heat Rate", title = "Nov 29, 2018 - Feb 27, 2019 (Furlough Year)") +
  scale_x_date(limits = c(as.Date("2018-11-29"), as.Date("2019-02-27")),date_labels = "%m-%d", breaks = tick_dates) +
  coord_cartesian(ylim = c(5, 50)) +
  geom_vline(xintercept = as.numeric(as.Date("2018-12-29")), linetype = "dashed", color = "red") +
  geom_vline(xintercept = as.numeric(as.Date("2019-01-27")), linetype = "dashed", color = "red") +
  theme(text = element_text(family = "Helvetica", size = 12)) +
  theme(panel.grid = element_blank(),
        axis.line = element_line(),
        panel.border = element_blank(),
        panel.background = element_rect(fill = "#F0F0F0"),
        axis.text = element_text(color = "#333333"),
        axis.title = element_text(face = "bold", size = 14, color = "#333333"),
        plot.title = element_text(face = "bold", size = 16, color = "#333333")) +
  theme_classic()
